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ABSTRACT 



This investigation develops an axi symmetric heat transfer- 
combustion model of a porous medium within a circular cylin- 
der. System flow is governed by Darcy's law. Carbon and 
air properties are treated as variables of temperature. A 
combined continuity-Darcy equation, an oxygen mass balance 
equation, and energy balance equations (one each) for air 
and carbon, describe the conservation laws of the system. 
Transport mechanisms for oxygen mass transfer are molecular 
diffusion and convective transport, and an oxygen consumption 
term to account for combustion is included. Heat transfer 
mechanisms included in the model are conduction and convec- 
tion. Radiation is accounted for at applicable boundaries 
only. Nonvolatile combustion is accounted for in the carbon 
energy and oxygen mass balance equations as a heat generation 
term of Arrhenius type. The numerical solution of four 
coupled, nonlinear, transient partial differential field 
equations is accomplished using the Galerkin formulation of 
the Finite Element Method. The effect of porosity on system 
behavior is examined. 
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INTRODUCTION 



I . 

A. PRIOR INVESTIGATIONS 

The problem of characterizing physical systems involving 
combustion and quantifying the accompanying thermal response 
has received attention in recent years. The level of diffi- 
culty encountered in a problem of this type precluded 
analytical as well as numerical solutions. The difficulty 
of the problem lies in the number of disciplines that encom- 
pass it. The problem involves the kinetics of combustion, 
heat and mass transfer mechanisms and fluid flow. Numerical 
solutions with increased efficiency in computation are now 
possible with high speed computers. 

The combustion problem has applications in the areas of 
forest fire control, energy conservation, underground 
(nuclear) fuel storage and waste disposal, and natural gas 
fire control. Major contributions may be made in the area 
of energy production and conservation by the analysis of 
heat generation and ignition characteristics of combustible 
materials . 

The characterization of combustion and heat transfer in 
porous media has been of particular interest. A brief survey 
of some of the investigations in the area of combustion and 
heat transfer in porous media follows. The intent is to 
acquaint the reader with work in the field that offered 
insight into the formulation of this investigation. 
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Investigations by Kordylewski on the influence of aero- 
dynamics on the critical parameters of thermal ignition 
[Ref. 1] , examined homogeneous combustion in a two-dimensional 
cylindrical reactor. The transient analysis involved 
Arrhenius combustion of a solid fuel. Heat transfer mechanisms 
considered were convection and conduction. Because the inves- 
tigation dealt with thermal ignition theory, reactant con- 
sumption was omitted. The flow field was assumed to be 
steady prior to ignition and constant fluid properties were 
assumed . 

Safety dictates the assessment of the structural integrity 
of a building after a severe fire. In order to predict 
stresses due to fire in buildings, the thermal response must 
be known. Sahota and Pagni [Ref. 2], formulated a transient 
solution for two-phase, two-component flow in one-dimensional 
porous concrete structures. The mechanisms considered in- 
cluded: heat conduction, molecular diffusion of gaseous 

components, and pressure-driven convective flow subject to 
Darcy's law. 

A sudden reduction in feed temperature in a packed-bed 
reactor leads to the transient temperature rise known as 
"wrong-way behavior." This behavior was investigated by 
Mehta, Sams, and Luss [Ref. 3]. The work identifies the 
important rate processes and parameters which cause the 
behavior, and generates an expression for predicting the 
magnitude of the maximum transient temperature. 
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Ignition parameters in porous solid fuels have been 
analyzed by Kim and Chung [Ref. 4], They investigated 
three porous solid fuel geometries (a semi-infinite slab, an 
infinitely long circular cylinder and a sphere) with con- 
stant energy and gaseous oxidant fluxes at the fuel surface. 
Laplace transformation of the nondimensionalized oxidant mass 
equation and fuel energy equation allowed for asymptotic 
solution of a nondimensionalized transient temperature 
equation which is valid in the neighborhood of the fuel 
surface. Observations included shorter ignition times for 
spheres compared to slabs. Times for ignition increased with 
fuel size and approached values of semi-infinite slabs 
asymptotically. Other effects of size and geometry of porous 
solid fuels on ignition parameters are presented. 

Saatdjian and Caltagirone [Ref. 5] investigated a transient 
two-dimensional combustion model with natural convection. 

A porous medium undergoing exothermic combustion was saturated 
by a gas and bounded by two impermeable planes. Permeability, 
fluid viscosity, thermal conductivity of the porous medium, 
and the heat of reaction were assumed to be constants. 

Horne and O'Sullivan [Ref. 6], Hickox [Ref. 7] and Chan 
and Banerjee [Ref. 8], investigated the natural convection 
phenomenon in porous media. Exothermic reactions were not 
addressed in these investigations. Horne and O'Sullivan 
considered the effects of variable viscosity on the stability 
of a porous layer in a transient two-dimensional problem. 
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The boundary configuration for this work was a two-dimensional 
region with a line of vertical symmetry, enclosed on both 
sides by impermeable, insulated sidewalls, and heated along 
half of the base. The convective flow of fluid through 
porous media heated from below has applications in the study 
of behavior of geothermal systems. Hickox studied convection 
as an application to sub-seabed disposal of nuclear waste. 

The problem was a two-dimensional transient analysis of free 
convection produced by a concentrated heat source (implanted 
container of waste material) in the subsurface sedimentary 
layer of a seabed. The porous medium was assumed to be rigid, 
homogeneous and isotropic. Density changes of the fluid were 
accounted for only in the buoyancy term of Darcy's law. 
Permeability, viscosity and thermal conductivity were assumed 
to be constant. Chan and Banerjee conducted a transient 
three-dimensional analysis of natural convection in porous 
media. In addition to convective heat transfer, conduction 
was also considered between solid spherical particles that 
comprise the porous medium. The porous medium was considered 
homogeneous and isotropic in its physical properties including 
permeability and thermal conductivity, both of which were 
assumed to be constant with temperature. Fluid density was 
considered to be constant except in the buoyancy term of 
Darcy’s law. 

In 1982, Vatikiotis [Ref. 9] considered a transient one- 
dimensional heat transfer and combustion model of a porous 



17 



medium. The heat transfer mechanisms included were conduc- 
tion, convection and radiation. One problem considered was 
a porous mat subject to Arrhenius combustion in which all 
properties were temperature dependent. The capillary serial 
(permeability) model introduced by Scheidegger [Ref. 10] was 
employed . 

The current investigation considers the effects of 
porosity on system behavior for a two-dimensional (axi- 
symmetric) model of combustion in a porous medium. A great 
part of the following analysis may be found in Vatikiotis 
[Ref. 9] . The pertinent parts are repeated here for conven- 
ience to the reader. Deviations are cited and are for the 
most part due to the axisymmetric two-dimensional aspect of 
the present model vice the previous one-dimensional model. 

The storage scheme employed by the numerical model warrants 
the reader's attention. The savings in computer storage 
realized by utilizing the method of Franke and Salinas 
[Ref. 29], is substantial and is addressed in the Numerical 
Section of this work. The nonlinear combustion term is 
treated as . a bilinear spatial operator of carbon temperature 
and oxygen concentration, in contrast to the Vatikiotis treat- 
ment of the term as an excitation vector. The idea behind 
the present treatment was to capture the effects of the non- 
linear combustion term on both of these dependent variables. 

It was felt that this would alleviate some numerical diffi- 
culties. A numerical model is formulated and results are 
presented. 
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II. PROBLEM DESCRIPTION 



The problem under investigation is that of combustion 
of a porous fuel (carbon) imbedded in a cylindrical container. 
Pressure gradients induce convective currents through the 
medium. The air flow produces two opposing effects: 

(1) internal convective heat transfer, and (2) a supply of 
oxygen to promote heat generation through combustion. 
Extinguishment or sustained combustion of the porous medium 
depend on the interaction of these effects. If heat transfer 
dominates, the combustion will proceed to extinguishment, 
otherwise combustion will continue. A mathematical model was 
developed to provide an understanding of this interaction 
and its effect on thermal behavior. 

The mathematical model is formulated as follows. Energy 
balances on the carbon and convected air provide heat 
transfer equations for each. The heat transfer mechanisms 
incorporated in the model are: (1) conduction, (2) convection 

and (3) radiation (where applicable, at boundaries) . In 
addition, nonvolatile combustion is included in the carbon 
heat transfer equation as a fractional order heat generation 
term of Arrhenius type. 

The conservation of species law is applied to the oxygen 
molecule concentration to yield a third equation. The 
resulting oxygen mass transfer equation includes the 
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transfer mechanism of: (1) molecular diffusion, and 

(2) convective transport of species. An oxygen consumption 

term due to combustion is included. 

The fourth field equation involves the system pressure 
gradient. The equation is a combined Darcy's law and air 
mass continuity equation for the system. 

The conservation equations describing the system field 
are four transient, coupled, nonlinear partial differential 
equations. The four equations are solved by a two-dimensional 
Galerkin formulation of the Finite Element Method. The inte- 
gration scheme used for this highly stiff system is one 
presented by Gear and modified by Franke [Ref. 30]. The 
scheme is especially suited to systems of implicit and stiff 
differential equations. The integration scheme is used in 
conjunction with an optimal compact storage scheme proposed 
by Franke and Salinas [Ref. 29] . 
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III. 



THEORY AND BACKGROUND 



A. DESCRIPTION OF THE POROUS MEDIUM 

In this work, a porous medium is considered to be a 
solid containing interconnecting pores that allow fluid to 
permeate and flow through the solid. Either of two classes 
of porous media, consolidated (solid and rigid) and uncon- 
solidated (comprised of discrete particles as found in 
granular beds) may be considered. Each class of porous 
media may have isotropic nonhomogeneous properties. 

The common characteristics of all porous media are: 

(1) porosity, (2) specific internal area, (3) pore diameter, 
and (4) tortuosity. Porosity, p, is defined as the ratio of 
void volume to total volume. The specific internal area, Z, 
is the ratio of internal surface area to bulk volume. In 
general, the distribution of pore size in a porous medium is 
random (nonhomogeneous) and dynamic (subject to small 
strains induced by the transient pressure field) . The situa- 
tion motivates the investigator to treat the porous medium 
as a continuum possessing idealistic geometrical properties 
of porosity and pore diameter. The specific internal area 
is obtained from the porosity model. Though many conventions 
exist for the definition of a conceptual pore diameter, in 
this study a hydraulic diameter analogy i.s employed. The 
tortuosity, x , of a porous medium is the ratio of the flow 
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path to the straight (line) path displacement of the 
particle. Permeability, a measure of hydraulic conductivity, 
is a property of the porous medium that depends upon the 
four characteristics mentioned above. Scheidegger [Ref. 10] 
presents methods used to measure the properties of porous 
media. Methods discussed are essentially experimental in 
nature . 

The porous medium was modelled as shown in Figure 3.1. 

In Figure 3.1, D is the particle center- to-center distance, 
and d is the particle diameter. From the idealized geometry, 
the porosity for spherical particles is. 



P 



- -(-) 3 
6 ( D 



(3.1) 



The pore diameter is obtained from an expression proposed 
by Carman [Ref. 11], 



6 




(3.2) 



Equation 3.2 is analogous to a more familiar form of mean 
hydraulic diameter, 4v/A, where v is the void volume and A 
is the wetted surface area. The specific internal area, Z, 
based on the idealized geometry of Figure 3.1 may be expressed 
as , 



Z 



2 V 



(3.3) 
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for spherical particles. Equation 3.3 is sometimes used 
and assumes one-half of the total internal surface area 
is effective for convective heat transfer. The fractional 
amount of total area is an estimate based on Fontenot's [Ref. 37] 
experimental results and does not generally apply to porous 
media [Ref. 9]. The Kozeny relations, alternate expressions 
of the specific internal area, are discussed by Scheidegger 
[Ref. 10] . Advantages of the Kozeny relations are fair 
agreement with experimental values and calculations that are 
independent of particle shape. A disadvantage is the 
failure of the relations to predict accurate values of z at 
high values of porosity [Ref. 9]. For the geometric 
configuration of Figure 3.1, the tortuosity depends on the 
ratio d/D. Carman [Ref. 11] presents a table of measured 
tortuosity factors of various materials and geometries, and 
points out differences between analytical determinations of 
tortuosity. In this study, his recommended value of 1.4 is 
used. Particle size decreases with consumption. As a 
result, thermophysical properties which depend on particle 
diameter, as well as temperature, are functions of time- and 
space. The numerical model presented assumes matrix 
rigidity as particle diameter decreases. This assumption 
gives rise to an increase in porosity during combustion. 

Although Scheidegger [Ref. 10] in his discussion of the 
packing theory of spheres, reports .875 porosity as the 
threshold for stability in a porous matrix. Carman [Ref. 11] 
reports on investigations performed on porosities as high 
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as 0.99. Changes in particle diameter are incorporated in 
the model and are discussed in Section III.D with carbon 
combustion. 

B. DARCY'S LAW AND PORE VELOCITY 

The Reynolds number for porous media is defined by 

p s d 

Re = — (3.4) 

U 

where s is the local pore velocity, p is the mass density of 
air, and y is the dynamic viscosity. The magnitude of the 
Reynolds number indicates whether fluid motion is dominated 
by molecular, viscous, or inertial effects. Most investi- 
gations of flow through porous media indicate flow regimes 
of viscous and inertially dominated flows. The Navier- 
Stokes equations apply to fluid motion possessing such 
Reynolds numbers. Because of the geometry involved in a 
consolidated (rigid) porous medium and the no-slip boundary 
condition (i.e., s = 0 at a solid-fluid interface), solution 
of the Navier-Stokes equations is difficult for a porous 
medium. Scheidegger [Ref. 10] points out that extensive 
experimental work with porous media indicates fluid flow is 
governed by Darcy's law for the range of Reynolds numbers 
where viscous effects dominate. The upper limit (velocity) 
Reynolds number in these investigations is subject to dis- 
agreement and varies from 0.1 to 75. Inertial effects 
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increase with increasing Reynolds numbers. Boffa [Ref. 12] 
has shown that for a fixed Reynolds number, inertial effects 
diminish with increasing air temperature. Darcy's law for 
two-dimensional flow is, 



Q 



-m , dP . / dP 

— • (3— r +(j— 
y dr dz 



- P 



P 3 g) z) 

a. 



(3.5) 



where Q is the filter velocity or the volumetric flow rate 
per unit cross-sectional area; m is the specific permeability; 

/v /\ 

g is the gravitational acceleration; r and z are the 
unit vectors in the r and z directions, respectively; and 
dP/dr and dP/dz are the pressure gradients in the r and z 
directions, respectively. The assumption here is that the r 
(radial) and z (axial) velocity components each react to the 
pressure independent of the other. The specific permeability 
of the porous medium used in the model is. 



m 



_2_(1 

96 t 



6, 2 



(3.6) 



Expression 3.6 is based on a cappilaric-serial model given 
by Scheidegger [Ref. 10]. Expressions for permeability vary 
with the physical assumptions of the flow paths incorporated 
into the model. Permeability also varies with the pore size 
distribution assumed. Physically, permeability and porosity 
are not related [Ref. 9]. Porosity is a quantifiable 
property of the porous medium, whereas permeability is a 
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a property specified by a 



constitutive property (i.e., 
constitutive relation, Darcy's law). The specific permea- 
bility is proportional to the filter velocity induced by a 
unit pressure gradient. Darcy's law is not derived from 
first principles but through exhaustive experimental 
analyses . 

The Dupuit-Forcheimer assumption, addressed in Carman 
[Ref. 11], relates the local pore velocity to the filter 
velocity by. 



The hypothesis is that the local pore velocity is greater 
than the filter velocity. The actual velocity in a single 
pore is a function of the fluid element's location within 
the pore. The Dupuit-Forcheimer assumption defines an 
"average" velocity within the pore. 

The continuity equation for two-dimensional flow in 
porous media with nonconstant porosity distributions is. 



Invoking the Dupuit-Forcheimer assumption and Darcy's law, 
the continuity equation becomes, 



Q 



P V 
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Equation 3.9 (neglecting body forces) is one of four field equations 
cast into a finite element formulation later in this work. From the 
pressure distribution, the pore velocity distribution is 
obtained by invoking Darcy's law and the Dupuit-Forcheimer 
assumption. The derivation of Equation 3.9 is presented in 
Appendix A. 

C. SEMENOV MODEL OF COMBUSTION 

The model of Semenov described in Frank-Kamenetskii 
[Ref. 13] and Vulis [Ref. 14], is the combustion model 
employed herein. The basis of the model is the relation of 

t 

reaction rate to temperature and the interaction of heat 
generation and heat transfer. The reaction rate expression 
R , is the Arrhenius expression for a simple n-th order 
reaction. 



* a n , -E x 
= A <{> exp ( — — ) 

R T 
u c 



(3.10) 



where A is the time constant of the chemical reaction, E is 
the activation energy, R^ is the Universal Gas Constant, 

A 

T c is the absolute carbon temperature, and 4> is the oxygen 
concentration. In a simple reaction, the reaction rate 
depends on the concentrations of reactants and not on the 
products. The heat generated by the reaction is obtained 
by multiplying Expression 3.10 by the heat (enthalpy) of 
combustion. As explained in (Ref. 13], a theory of combustion 
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can be constructed only if certain reasonable assumptions 
are made. For example, in order to regard the initial 
states as stationary, one must neglect the reaction rates 
at these states. The empirical law of chemical kinetics 
embodied in the Arrhenius expression tells us that the reac- 
tion rate never goes to zero but falls off exponentially 
with a decrease in temperature. Neglect of the reaction 
rate at the initial states is necessary to achieve initial 
equilibrium. In a finite range of temperatures above the 
initial states, neglect of reaction rates is also justifiable 

[Ref. 13]. At room temperature the reaction rates are on 

2 

the order of l.E-15 lbm-carbon/f t -hr or smaller, vice 

2 

l.E + 3 lbm-carbon/f t -hr at 1500 °F. 

D. ARRHENIUS LAW OF REACTION RATES 

In 1934, Parker and Hottel [Ref. 15] proposed an 
Arrhenius expression for the reaction rate of carbon reacting 
in air. The expression assumed a simple first order reaction 
for the combustion of carbon in air. Frank-Kamenetskii 
[Ref. 13] has shown that the Parker and Hottel data is better 
correlated by a fractional order reaction, n between 1/3 
and 2/3. A reaction order of 1/2 yields, 

R = 2.065 x 10 6 (R n (J>) 1//2 exp(-57,240/R -T ) (3.11) 

C U 2 u c 

2 

In Expression 3.11, R c is in units of lbm-carbon/f t -hr, 

R q is the gas constant for oxygen (48.29 f t-lbf/lbm-R) , 
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4> is in lbm/ft^, R u is 1.986 Btu/lbmole-R, and T c is in 
Rankine. In order to determine the rate of heat generation 
and the rate of oxygen consumption, the chemical reaction 
for the combustion process must be considered. Although many 
complex chains in chemical kinetics generally describe 
combustion, a simple two-product analysis is employed in this 
formulation. For nonvolatile combustion of carbon and 
oxygen, two reactions that describe the process are, 



C + 2 ° 2 



CO 



(3.12) 



C + o 2 - co 2 



(3.13) 



where 0 denotes oxygen; and C, CO and C0 2 denote carbon, 
carbon monoxide and carbon dioxide, respectively. The ratio 
of the mass rates of carbon monoxide to carbon dioxide 
produced increases with increasing temperature. Arthur 
[Ref. 16] presents an expression for the rate ratio as a 
function of temperature (in Kelvin) . 




2500 exp(-6240/T c ) 



(3.14) 



The expression is valid for temperatures between 790 and 
1170 K (1310-2110 degrees Fahrenheit). As a result of this 
temperature dependency, the stochiometric ratio and the heat 
of reaction are functions of temperature. Defining the 
fraction of carbon monoxide being produced by/ 
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(3.15) 



? = (-££-) / (1 + (-££-) ) 
CO C0 2 C0 2 



and the fraction of carbon dioxide as. 



CO, 



1/(1 + 



the heat of combustion is then expressed as, 



(3.16) 



AH 



R 



F co AH co + F co. 



AH 



CO, 



(3.17) 



Values for the heats of combustion, AH and AH , as 
functions of temperature may be obtained from JANAF (Joint 
Army, Navy, and Air Force) tables [Ref. 17]. For the range 
of temperatures being investigated (80-2000 degrees Fahren- 
heit) , the heats of combustion are 3966.3 Btu/lbm for F 

and 14,121. Btu/lbm for F . Frank-Kamenetskii [Ref. 13] 

C°2 

points out that over narrow ranges of temperature it is 
permissible to use an approximate expression which correctly 
describes the reaction rate. The stochiometric ratio (fuel 
to oxygen) of the overall reaction is. 



■R 



f co f co. 



/(f 



CO- CO 



f F 
CO CO, 



(3.18) 



where f CQ is the stochiometric ratio for the reaction 3.12 

and f is the stochiometric ratio for the reaction 3.13. 

The rate of heat generation, R , and the rate of oxygen 

9T 
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consumption may now be expressed by, 



R 



g 




(3.19) 



R 



0 



2 




(3.20) 



Parker's and Hottel's work [Ref. 15] and Arthur's work 
[Ref. 16] were conducted with specific types of carbon. 

Tables and references exist in Smoot and Pratt [Ref. 18] and 
Frank-Kamenetskii [Ref. 13] for rate expressions using other 
types of carbon. The present model allows for any simple 
fractional order rate expression of Arrhenius type to 
account for carbon consumption with CO and CO 2 byproducts. 

For this work the Parker and Hottel rate expression as modi- 
fied by Frank-Kamenetskii (n = 1/2) is used. 

Particle diameter decreases with progressive combustion. 
The rate of decrease depends on the amount of carbon con- 
sumed at a point over time. Observations have shown that 
the effect of decreasing particle diameter is significant 
when the reaction is concentrated in a small region of the 
porous medium. To account for this, an expression for the 
time rate of diameter change as a function of reaction rate 
is derived (Appendix B) . For spherical particles the equation 
is. 



O 



-2R Z D 3 / ( iTp d 2 ) 

V-* 



d 



(3.21) 
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where p c is the bulk mass density of carbon. The diameter 
and the reaction rate are functions of time and space. 

E. CARBON HEAT TRANSFER EQUATION FOR POROUS MEDIA 

The heat transfer equations of the present investigation incor- 
porate: (1) radiation (at boundaries where applicable), (2) internal 

convection, (3) conduction, (4) internal combustion, (5) temperature 
dependency of properties, and (6) compressibility effects 
of air into a two-dimensional (cylindrical coordinate) 
formulation. The carbon energy conservation equation is, 

3T 

(l-P)P c C c = V- (1-p) (k e ) (VTJ - hZ(T c -T a ) + R g Z (3.22) 



The derivation of Equation 3.22 is presented in Appendix A. 
The effective conductivity, k g , of the porous solid was 
proposed by Russel [Ref. 24], 



k = k 
e a 



P ,2/3 + j^-(l - P ,2/3 ) 

c 

,2/3 , ^ k a , . ,2/3 ^ , 

P - P + — (1 " P + P 

c 



(3.23) 



where k and k are bulk thermal conductivities of carbon 
c a 

and air respectively, and p' = 1-p. Russel's expression, 
which is based on an electrical analogy, is valid for the 
full range of porosities, 0.0 to 1.0. 

Because of the difficulties encountered in a radiative 
analogy to Fourier's law of conduction, the particle to 
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particle radiative exchange in the porous medium is omitted. 

The difficulties are: (1) the geometry does not easily allow 

one to derive an expression for a "sink" temperature to use in 
a linearized approximation of the Stef an-Boltzman equation, and 
(2) the multi-wavelength characteristics of the radiation 
phenomenon are not easily incorporated. 

F. HEAT TRANSFER EQUATION FOR AIR IN POROUS MEDIA 

The internal convective heat transfer coefficient of 
Yoshida, Ramaswami , and Hougen [Ref. 33], is given by, 

h = 0.91 Re ' ' 51 C G(C y/k ) _2/3 ]f (3.24) 

a a a m 

where ip is equal to 1 for spherical particles, G is a pseudo 

mass velocity given by p p s and C is the specific heat 

a a 

of air at constant pressure. The f subscript indicates 
properties are to be evaluated at film temperature. Re' is 
a pseudo Reynolds number defined by, 

Re' = G/z y ip (3.25) 

The air properties, as well as the internal convection coeffi- 
cient, h, are temperature dependent. The reaction rate in 
Equation 3.22 is given by Expression 3.19. 

An energy balance on the air within the porous medium 
obtains the second heat transfer equation, 
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(3.26) 



3T 



PP C 



a a 3 1 



V- (p k VT ) + hZ (T - T ) 

3 a. C 3 



- p p C (V* V) T - (V-V)p P 

3 a. a. 



The density of air, p , is approximated by the ideal gas law. 

cl 

Pressure terms are due to the compressibility of air. The 
derivation of Equation 3.26 is presented in Appendix A. All 
properties in Equation 3.26 are temperature dependent. The 
properties of standard air were used in the model. Vatikiotis 
[Ref. 9] points out that tolerable differences (average of 
7% difference) between standard air properties and properties 
accounting for the presence of byproducts CO and CO^ / is 
acceptable. Increased accuracy would introduce an additional 
mass balance equation for either CO or CC^ • Polynomial ex- 
pressions used to calculate the thermophysical properties of 
air are those presented by Vatikiotis [Ref. 9]. The expressions 
are presented in Appendix B. 

G. OXYGEN DIFFUSION EQUATION FOR POROUS MEDIA 

The fourth field equation necessary to complete the 
system of equations is provided by an oxygen species conser- 
vation requirement. Transport mechanisms included in the 
model are molecular diffusion (Fick's law), convective mass 
flow and oxygen consumption due to combustion. Pressure and 
temperature induced concentration gradients are considered 
negligible. Vatikiotis [Ref. 9] , provides examples for which 
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diffusion arising from pressure and temperature gradients 
is important. The oxygen mass balance equation is, 

P |f = V- (p T? e V<J>) - V- (p <{>V) - R Z (3.27) 

The derivation of Equation 3.27 is presented in Appendix A. 

The effective diffusivity, V , proposed by Denbigh and 
Turner [Ref. 31] for a porous medium is 

V q = P/t (3.28) 

Expression 3.28 accounts for the tortuosity encountered by 
the oxygen molecules as they flow through the porous medium. 

The semi-empirical expression proopsed by Gilliland [Ref. 34] 
is used to obtain the diffusivity of oxygen into air. The 
expression is, 

V = 435.7 T^ /2 (M~ 1 +m" 1 ) 1/2 /[P(V^ /3 +V^ /3 ) 2 ] (3.29) 

3 a. ^ 2 3 ^2 

2 

where V is m units of cm /s, P is the total pressure m Pa, 

V and V are the molecular volumes of air and oxygen, 

3. U 2 

respectively. M and M are the molecular weights of air 

a u 2 

and oxygen. The values of V and V may be obtained in 
Holman [Ref. 35] as 29.9 and 7 . 4 cm , respectively. The 
oxygen consumption term in Equation 3.27 is given by Expression 
3.20. 
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H. BOUNDARY CONDITIONS 



Boundary conditions employed were as follows for the 
carbon , 



3T 

c 

3r 



— = 0 
r 

o 



= 0 



C 



u-p) (k e ) || 



— = 0 
z 

o 



= -q 



c 



where q g is the starting heat flux. 



^ rp 

(1-P) (k e ) gj 



— = 1 



~4 ^ 4 

-OE (T - T ) 

c CO 



c 



(1-p) (k ) 

^ e 3r 



-£- = i 
r 

o 
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C 



^4 "4 

h (T -T ) + os (T -T ) 
r c oo c =o 



For air. 



(2-D problems) 
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3T 
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(3.35a) 



T 

a 




0 



= T 



OO 



P k 



3T 
a 

a 3z 



= 1 



- PP C v (T - T ) 

a a. a. Q o 



(3.36) 



3T 



P k 



a 3r 



r 



(1-D problems) 



(3.37) 



. , ~4 ~4 

h r (T a " T J + ae (T " T ) 

X d °o 3 oo 



(2-D problems) 



For pressure. 



3p 

3r 



— = 0 
r 

o 



= 0 



= P. 



= 0 



= P, 



= 1 



3p 

3r 



o 



= 0 



(3.38) 



(3.39) 



(3.40) 



(3.41) 
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Finally, for 0 2 concentration, 



34 > 

3r 




0 



0 



(3.42) 



P V 



e 




v(<p - (p^) 



(3.43) 



4 > 



z 

z 

o 



= 0 



= <f> 



(3.43a) 



P V 



e 




1 



v ((f) - <f> ) 
00 



(3.44) 



P V 



e 



3(f) 

3r 




(3.45) 



Equations 3.35, 3.36, 3.43, and 3.44 correspond to convective 

flux conditions on air temperature and 0 2 concentration. 

At the air inlet (— = 0.0), Dirichlet conditions. Equations 

o 

3.35a and 3.43a may be imposed on the air temperature and 
oxygen concentration. The idea behind this treatment is that 
in the presence of a semi-infinite medium (ambient air) , the 
air and C> 2 concentration may be considered, to a first 
approximation, very near ambient conditions. Equations 3.30, 
3.34, 3.38 and 3.42 correspond to symmetry conditions at 
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— = 0.0. Equations 3.33, 3.37, 3.41 and 3.45 correspond 
o 

to impermeable boundaries and insulated boundaries at = 1.0 

and are used in conjunction with one-dimensional problems. 

The second part of Equations 3.33 and 3.37 are the boundary 

conditions that represent a relaxation of the radial insulation 

of carbon and air temperatures which allow the system to 

become a two-dimensional heat transfer problem. It must be 

pointed out however that the heat transfer coefficients of 

Equations 3.33 and 3.37 are not easily obtained. This work 

could not proceed if the following assumptions were not made. 

It was assumed that the air and carbon had near equal profiles 

at — = 1.0 so that the same heat transfer coefficient would apply 
o 

to both variables. If Equation 3.35 is used as a boundary 
condition, the air temperature follows closely the carbon tempera- 
ture at — = 1.0. Furthermore it was assumed for simplicity of 
r o 

calculation that the heat transfer coefficient was constant. 

The difficulty is as follows. In order to determine the 

heat transfer coefficient for the impressed temperature 

profiles at — = 1.0, the temperatures themselves must be 
o 

known. The correct procedure would involve an initial esti- 
mate of the profile, calculation of the heat transfer coeffi- 
cient, solution of the problem for one integration and 
verification of the assumed and calculated temperature 
profiles. Upon suitable verification of the profiles, the 
same procedure is performed for the next integration. Other 
problems arise for example if the radial heat transfer at 
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— = 1.0 is of natural convective type. At the exit — = 1.0, 
o z o 

functional and derivative continuity must be demonstrated 

for pressure, oxygen concentration and air temperature at 
the vertical boundary separating the stream of exiting air 
and the convective boundary layer. This is known as a conju- 
gate problem. The above method is a realistic method for 
the treatment of this problem but admittedly it is beyond the 
scope of this work. The numerical code has provisions for 
incorporating an average of isoflux and isothermal natural 
convective Nusselt number formulations (Churchill correlations) 
for vertical cylinders obtained from [Ref. 36]. Although the 
subroutine has been written there has been no attempt to 
date to implement it. It is known that physically heat 
transfer coefficients lie between those generated from iso- 
thermal and isoflux considerations. However, simply having 
empirical correlations for evaluating the heat transfer 
coefficients does not eliminate the need of the iterative 
scheme described above. Additional considerations must be 
addressed in this context. In order to effect an impermeable 

vertical boundary at = 1.0, the system must be enveloped 

o 

by a cylindrical container fabricated of some type of material. 
This material has an inherent potential to absorb internal 
energy that transits from the initial system to the boundary 
where some flowing medium is able to convect energy away 
(Figure 3.2). The material's ability to absorb energy adds 
a thermal resistance to the heat transmission in an electical 
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analog sense (Figure 3.3) and thus affects the temperature 
profile and ultimately the heat transfer coefficients. Of 
necessity then, the material's inherent ability to absorb 
energy is denied. No attempt is made here to conjecture as 
to how physically one may impose the radial boundary conditions 
implemented below. If it were possible to rigorously apply 
the heat transfer boundary condition described above, the 
numerical model would be capable of obtaining a solution. 

The model was tested at various values of constant heat trans- 
fer coefficient: 0, 1, 5, 50, 100, 500, to determine the 

effect, of the heat transfer coefficient on the solution. It 
was found that until the heat transfer coefficient gets 
significantly large compared to the start flux applied to 

the carbon at the air inlet boundary (in this case above 
2 

100 Btu/ft -hr compared to 1500 flux applied to carbon) , the 

two-dimensional effects on the temperature profile were 

localized in the — = 0.85 to 1.0 region. The temperature 

o 

profiles for — less than 0.85 were essentially one-dimensional 
o 

(constant value independent of r) . 

I. INITIAL CONDITIONS 

The model is developed so each problem begins at a 
uniform initial temperature. A heat flux applied to the 
carbon along a boundary provides a means of bringing the 
system to a temperature level where reaction terms are 
sufficiently high to generate combustion. The heat flux is 
treated as follows, 
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(3.46) 



(l-p)?(k ) 
e 





-q 



it 

s 



Expression 3.46 is incorporated into the numerical formula- 
tion as Neumann boundary condition for carbon temperature. 
The initial heat flux may be turned off at any specified 
time. The boundary condition of Expression 3.46 then 
switches to a Cauchy (mixed) boundary condition to account 
for radiative heat transfer from the carbon particles at the 
boundary to the ambient air. Convective heat transfer from 
the carbon particles to the air is treated through the 
internal heat transfer coefficient. The above procedures 
for treating problem initiation obviate the need of trying 
to specify initial conditions which may in general be 
arbitrary for each problem. 
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Figure 3.1 Geometric Model of Porous Medium 
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Figure 3.2 



Top View of Radial Convection Circuit 
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Figure 3.3 Electrical Analogy of Thermal Circuit 
at r/r Q = 1.0 
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IV. NUMERICAL CONSIDERATIONS 



A. TREATMENT OF NONLINEAR COMBUSTION TERM 

The treatment of the nonlinear combustion term is now 
discussed in brief detail highlighting the advantages of 
each treatment. A more comprehensive discussion is presented 
in Appendix C. The interested reader is encouraged to 
review the detailed analysis. There are several ways to 
treat the Arrhenius reaction rate expression. 

1 . Excitation (Force) Term Treatment 

Perhaps the simplest method for incorporating the 
combustion term is to evaluate it at each time step and use 
the evaluated value (held constant) as an estimate of the 
mean value during the next integration. It is realized that 
in fact, the value is not constant in real time. This treat- 
ment, however, provides a means of incorporating the term 
into the system of equations as a first approximation with 
relatively little computational effort. Averaging techniques 
exist for improving upon this method. In this scheme the 
excitation vector is modified at each nodal point in the 
carbon temperature and oxygen concentration equations to 
include the combustion terms. The exponential reaction rate 
term that appears in both the porous solid and oxygen diffu- 
sion equation is, 



A 

R c = A <p n exp (-E/RT ) 



( 4 . 1 ) 
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In the carbon equation. Expression 4.1 appears as. 



R g = AH r • R c (4.2) 

In the oxygen concentration equation. Expression 4.1 appears 
as , 



R 



0 



2 




(4.3) 



The term that is to be evaluated at the last time step and 
to be held constant over the next time step is. 



R 



= {A (p 



n 



exp(-E/RTj } 



'n- 1 



(4.4) 



where the superscript indicates evaluation occurring 

at the previous time step. 

2 . Linear Operator Treatment of 0 ^ Concentration 

In order to realize an improvement over the first 
treatment one may retain a portion of the reaction expression 
as an operator by making the following rearrangement: 

/\ * t 

R c = {A <J) n-1 exp(-E/RT c ) } n_1 ♦ [<j>] (4.5) 



where [<j>] denotes a spatial operator treatment of the 
response variable. 
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3 . Bilinear Operator Treatment of 0 ^ and 



The bilinear operator treatment of the reaction rate 
is the present method of treatment. The reaction rate term 
is rearranged as follows, 



1 ^ 

A <p n X exp(-E/RT ) *t 1 
= { — } n 1 • [T <j>] 



(4.6) 



Letting 



T 

Hhi 



i = 1,3 



(4.7) 



and 



4> * ^e 4i , i = 1,3 (4.8) 

and invoking natural coordinates [Ref. 28], an elementally 
averaged contribution results in the following area integral, 

C = C / / 5 ? T 0 1 £ T 6 dA (4.9) 

e 

The double subscript permutation of carbon temperature and 
C >2 concentration results in a 3x9 elemental matrix that is 
distributed into the system equations. In this scheme at 
each nodal point within an element, carbon temperature and 
C >2 concentration equations receive nine terms arising from 
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combustion occurring within the element. An element loop is 
performed to account for combustion terms in an elementally 
averaged sense over the entire system spatial domain. 

B. OPTIMAL COMPACT STORAGE (OCS) 

The optimal compact storage scheme, presented by Franke 
and Salinas [Ref. 29], employed in the solution of the 
combustion problem, allowed the computational effort to 
proceed with substantial savings in computer storage, 
computer-run time, and ultimately dollar cost per run. In 
the optimal compact storage scheme, only the non-zero coeffi 
cients are stored in a vector array. This results in a 
very significant reduction of the storage area compared to 
banded storage. One might encounter problems on the order 
of 1000 x 1000 DOF. (In this problem there are four degrees 
of freedom at each nodal point. A 17 x 17 grid requires a 
1156 x 1156 matrix if full storage is employed.) For banded 
storage the bandwidth might be 200. The storage ratio for 
bandwidth to full storage would be 0.2. Using OCS, a con- 
servative estimate of the average number of equation entries 
in this model is 20. The ratio of OCS to banded storage 
is 0.1. Thus in this example, the savings realized by OCS 
vice full storage is 98% 1 

C. GRID CONVERGENCE 

A convergence study was done on the response field for 
several grids. It was found that non-uniform grids were far 
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superior to uniform grids. In fact, for many problems, 
computer storage limitations require that a non-uniform grid 
be used to obtain reasonably accurate solutions. Non- 
uniform grids allow the investigator to take advantage of 
the knowledge of areas of severe combustion induced gradients. 
"Stacking" elements in these areas assists the algorithm in 
producing more accurate and numerically stable results. For 
similar degrees of accuracy with uniform grids, it is esti- 
mated that grids on the order of 1000 * 1000 nodal points 
(and larger) would have been required. Nondimensionalized 
length discretizations on the order of .001-. 005 and smaller 
near the air inlet boundary permit excitations on the system 
boundaries to be propagated accurately through the medium. 
Until discretizations on the order of non-dimensionalized 
lengths of 0.001 were employed, numerical instability was 
encountered in the carbon temperature field. This instability 
was manifested by severe overshoot, leading to negative 
temperatures in close proximity to severe thermal gradients. 

For non-uniform grids it was found that grids below 10 
nodal points in each direction were not adequate. The 
temperatures obtained from these models were lower due to the 
relatively low degree of grid refinement or discretization 
possible with such few points. Grid refinement in this type 
of problem is essential to capture the high activity (large 
gradients) that is inherent in the combustion phenomenon. 

A convergence study was done on several nodal points at 
various times for non-uniform grids at 12 x 12, 17 x 17, and 
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22 x 22 nodal points. The solution changed 5% from a 12 x 12 
to a 22 x 22 nodal point model. Next, a 17 x 17 point non- 
uniform grid model was investigated. It was found that 
the solution changed approximately 2% from the 22 x 22 nodal 
point model. At this point, it was determined that the 
12 x 12 non-uniform nodal point afforded the desired balance 
between cost, computational effort and computer-run time. 





grid size 


kmax 


clock runtime (min. ) 




(n x n) 








12 


6,398 


70 




17 


13,223 


100 




22 


22,498 


110 


D. 


MODEL VALIDATION TESTS 






In order to validate the capabilities and accuracy of 


the 


numerical code. 


several tests 


were conducted. 




1. The Steady 


State Problem 





First, it was necessary to ensure the model's ability 

to recognize a steady state condition. This was the simplest 

of all tests (and first to be) performed. Initial (ambient) 

conditions were imposed uniformly throughout the system on 

the four fields: carbon temperature and air temperature 

(80 degrees F.), pressure (2,116.8 psf) and oxygen concen- 

3 

tration (0.0172 lbm/ft ). An initial integration time step 
equal to the minimum time step allowable (user input) l.E-6 hr., 
was selected. In less than five integrations the algorithm 
switched to the maximum allowable time step (user input) 
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5.E-1 hr., and iterated about the given initial fields. 

This sufficiently demonstrated the first of several tests 
requisite to assure numerical code compatibility with the 
integration scheme. 

2 . The Finite Slab 

Consider the slab of thickness 2L in Figure 4.1(a). 
One seeks the solution of the one-dimensional conduction 
equation , 



9T 

3x 



= a 



3 2 T 

^7 



subject to an initial condition, 



(4.10) 



T (x , 0 ) = T i (4.11) 

and uniform convection conditions at both surfaces: At 

x = ±L : 




± h (T - T ) 
o 



(4.12) 



The exact solution is given as a Fourier series in [Ref. 23], 




l 

i=l 



C . 
1 
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cos ( B i jt) 



(4.13) 



where , 



(4.14) 



4 sin 8^ 

C i 2 8^ + sin (28^) 

The constants 8^ are the roots of the transcendental alge- 
braic equation. 



3 • tan 8 • = B . = h L/k 

1 1 1 o 



(4.15) 



Thus the solutions have the Biot number of the slab as a 
parameter. The roots of Equation 4.15 are tabulated in 
Appendix A to [Ref. 25] . 

3 . The Cylinder 

Consider the cylinder of Figure 4.1(b). The suddenly 
immersed cylinder satisfies the equation/ 



3T _ a 3, 3T 
3t r 3r r 3r 



(4.16) 



subject to an initial condition. 



T (r , 0) = T ± 



(4.17) 



and a convection condition at the surface: At r = r , 

o 



-k 



3T 

3r 



h (T - T ) 
o o 



(4.18) 



The solution is given as a Fourier series in [Ref. 23] : 
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(4.19) 
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where , 



C. 

i 



w 



6 i Jq( 6 ;L ) + jJ(6 i ) 



(4.20) 



The constants 8^ are the roots of the algebriac equation, 



8 i Ji(8 i )/J 0 (8 i ) 



B . 
1 



h r /k 
o o 



(4.21) 



The functions J and J. are the Bessel functions of the first 

o 1 

kind whose numerical values are tabulated in most advanced 

engineering mathematics texts. For a limited range of 

interest, numerical values of J and J n are tabulated in 

o 1 

Appendix A of [Ref. 25]. The roots of Equation 4.21 are 
tabulated in [Ref. 23] . 

4 . Multidimensional Solutions by the Product Method 

The classic solutions for semi-infinite- and finite- 
thickness bodies (slabs, cylinders and spheres) may be used 
to generate solutions for multidimensional bodies. 

Use of the product solution technique is illustrated 
by Figure 4.2. The problem chosen to validate the heat 
transfer portion of the current model is the "sudden immer- 
sion" problem for a finite-length cylinder subjected to 
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uniform convection conditions (h ,T ) on all sides and 

o 00 

initial uniform temperature T^ . The differential equation 
to be solved is. 



1 a, do_. a 2 e = 1 _ae 

r 3r r 3r ^2 a 3t 



(4.22) 



where 0 = (T - T ) / (T^ - T ) . The initial condition is 



0 (r , x, 0) = 1 



(4.23) 



The convection boundary conditions are. 

At top and bottom: -k = ± h 0 (4.24) 

On the sides: -k 4-^- = h 0 (4.25) 

3r o 

The solution is the product of the two simpler analyses: 
the semi-infinite slab and the infinite cylinder. Let 

0 (x , r , t ) = e s iab^ X,1:) * 9 cyl^ r,t) (4.26) 

= P(x,t) • C (r , t) 

Substituting Equation 4.26 into Equations 4.22, 4.23 and 
4.24 and separating the variables, reduces the two-dimensional 
problem into two one-dimensional problems: 
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Slab : 



3 2 e 



3x 



, 36 

1 s 

a 3t 



(0 = 9 . , ) 

s slab 



6 (x , 0 ) = 1 



30 



-k 



3x 



x = ±L 



h 0 ( ±L , t) 

os ' 



1 a 96 

Cylinder: - ^(r -^) 



, 30 

1 c 

a 3t 



(9 =6 ,) 

c cyl 



0 c (r,O) = 1 



30 



-k 



3 r 



h 0 (r , t) 
o c o 



Thus Equation 4.26 is the exact solution to the sudden 
immersion problem of a right circular cylinder. 

5 . Validation Problem 

The validation problem is Example 4.9 in [Ref. 25]. 

The problem is stated here for the convenience of the reader. 

The short cylinder in Figure 4.3 is initially at 40° C 

and then plunged into a fluid with h = 300 W/m-K and T = 200° C, 

-9 2 

The material is bronze, k = 2 6 W/m-K and a = 8.6 • 10 m/s. 
Find the temperature at the center of the cylinder after 
5 minutes. 
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Solution: 176° C (348.8° F) . The details of the solution 

are found on pages 186-187 of [Ref. 25]. The product solution 

technique was utilized in obtaining a 5-term approximation 

2 

of the series solution for Fourier moduli (F q = at/L , 
equivalent to a nondimensionalized time) less than 0.2. For 
Fourier moduli greater than 0.2, a two-term approximation of 
the exact solution was used. Heisler [Ref. 26] points out 
that a one-term approximation of the exact series solution 
yields an accuracy to within 1% of the exact value for 
Fourier moduli greater than 0.2. Various values (5, 10, 20 
and 30 seconds) corresponding to values of Fourier moduli 
less than 0.2 were used in comparing model solutions with 
the exact (five-term approximations to the exact) solution. 
Values (3, 5, 7 and 9 minutes) corresponding to Fourier 
moduli greater than 0.2 were used to compare model solutions 
with the exact (one-terra approximations with 1% accuracy) 
solutions. Figures 4.4 through 4.8 show: (1) the effects 

of time on the solution of various grids, and (2) effects 
of grid refinement at equal values of time. Values plotted 
on the ordinate scales are deviation (percentage error) from 
the exact solutions versus time. Observation yields that for 
small values of time, grid refinement obtains more accurate 
results and obtains a faster convergence to steady state. 

6 . Testing the Model's Ability to Accept and Reject 
Heat 

An equivalent simple test in model validation is the 
model's ability to accept heat from a boundary flux and 
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reject it at a later time to its environment via convection. 

z 2 

The applied heat flux at — = 0 is 1500 Btu/ft -hr. 

o 

Figures 4.9-4.12 illustrate system response. Figure 4.13 
is the input data set used. 

7 . The One-Dimensional Problem 

Another step in the validation of the two-dimensional 
problem is an examination of a one-dimensional problem. One 
might infer that a two-dimensional model includes the one- 
dimensional model as a subset. This was demonstrated in a 
separate test. A one-dimensional problem may be imposed on 
the model in two ways. One way is to make the length to 
diameter ratio very large. (No restrictions on excitations 
is implied . ) 

Another method for eliciting one-dimensional behavior 

is to "excite the system in a one-dimensional fashion," 

z z 

i.e., apply boundary conditions at — = 0.0 and — = 1.0 

z z 

o o 

independent of radius and insulate radial boundaries at 

— = 0.0 and — = 1.0. The insulated boundary at — = 0.0 
o o o 

arises from an axisymmetric formulation of the problem. The 

zero gradient (completely insulated) boundary at — = 1.0 

o 

corresponds to a completely (all four variables) impermeable 

vertical boundary. For this problem, an adiabatic restriction 

is sufficient. In general, an adiabatic condition means an 

insulation of heat. In this problem, the implication is far 

greater because of the coupling between system fields. An 

adiabatic restriction at — = 1.0 prohibits non-zero radial 

r o 
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pressure gradients. Unconstrained radial pressure gradients 

at = 1.0 imply radial convection of air enthalpy and in a 
o 

non-isothermal fluid, this convection undermines the initial 

adiabatic assumption at — = 1.0. Oxygen gradients are also 

o 

negated by an adiabatic condition since oxygen convection 

by air must occur under the presence of nonzero pressure 

gradients at — = 1.0. 

o 

The two methods described above distinguish one-dimen- 
sional behavior arising from geometrical conditions and one- 
dimensional behavior resulting from restrictions on excitations. 
Figures 4.14-4.16 depict model one-dimensional solutions to 
"one dimensional excitations." Figure 4.17 is the input data 
set used. 
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Figure 4.1 Classical Transient 1-D Geometries 
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cylinder 
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Short cylinder 
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Product solution: 
eu, r. f) = P(x, t)C(r t t) 



Figure 4.2 Illustration of the Product Technique 
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Figure 4.3 Sample Problem, Cylinder Geometry 
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Figure 4.4 Grid Effects at Point A, (0,0) 
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Figure 4.5 Grid Effects at Point B, (0,1) 
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Figure 4.6 Grid Effects at Point C, (1,1) 
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Figure 4.7 Grid Effects at Point D, (1,0) 
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Figure 4 . 9 Model Accepting Heat (Carbon 
Temperature Profile) 
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Figure 4.10 Model Rejecting Heat (Carbon 
Temperature Profile) 
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Figure 4.11 Model Accepting Heat (Oxygen 
Concentration Profile) 
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Figure 4.12 Model Rejecting Heat (Oxygen 
Concentration Profile) 
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Figure 4.13 Accept/Reject Heat Problem Input Data Set 
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Figure 4.15 One Dimensional Problem (Carbon 
Temperature Profile) 
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Figure 4.17 One Dimensional Problem Data Set 
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V. DISCUSSION AND CONCLUSIONS 



Field problem behavior depends on three major factors: 

(1) Boundary conditions, (2) System excitation, and (3) The 
physical characteristics of the system. 

Time limitations did not permit this investigator to 
conduct an exhaustive analysis of each of these factors . 

Some preliminary analyses were performed to obtain some 
understanding of the effect of boundary conditions, and 
excitation on system behavior. 

There are a variety of physical parameters which govern 
system response, such as permeability, porosity, and the 
cylinder length- to-diameter ratio. Here only a brief 
investigation of the effect of porosity on system behavior 
was undertaken, and is reported. 

A. BOUNDARY CONDITIONS 

The boundary conditions used in the present investiga- 
tion were presented in Chapter III. Other boundary conditions 
are possible. For example, if in place of Eqs . 3.32 and 3.36, 

insulated boundary condition at — = 1.0 is imposed on the 

o 

carbon and air temperatures, then preliminary results indi- 
cate temperature response is higher for equal values of time. 
This behavior is due to the buildup and storage of energy 
within the system associated with an insulated boundary. 



an 
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Similarly, in place of Equation 3.44, a zero flux boundary 

condition at = 1.0 could have been imposed on the oxygen 
o 

concentration. The zero flux condition implies an impermeable 

boundary with respect to oxygen concentration fluxes. This 

boundary condition would apply to a very long cylinder (i.e., 

a regularity condition). In this investigation, a Cauchy 

(convective flux) boundary condition on the oxygen concentra- 

tion is imposed at — = 0.0 and 1.0. As oxygen is locally 

o 

consumed in the interior of the system, oxygen gradients are 
created. According to Fick's Law of Diffusion, a depleted 
oxygen region is replenished by the diffusion of oxygen from 
regions of relatively high concentration to regions of low 
concentration. Thus a convective flux boundary condition 
on oxygen causes the depletion of oxygen concentration to be 
retarded through the diffusion mechanism. 

B. EXCITATIONS 

The effects of high and low values of the heat flux 
2 

applied at — = 0.0, are predictable. High values of heat 
o 

flux lead to an accelerated development of both the carbon 

and air temperature profiles and 0 ^ concentration depletion 

2 

within the system. For high values of heat flux at — = 0.0, 

z o 

the numerical integration scheme eventually slows noticeably 
as a result of steep gradients observed at this boundary. 

C. PHYSICAL CHARACTERISTICS OF THE SYSTEM 

The physical characteristics of the system generate a 
significant effect in field problem solutions. In the 
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problems discussed here the transport properties k , k , 
m/y, V and the Arrhenius reaction rate term are governing 
factors in the respective field equations in which they 
appear. As previously discussed, the net balance between 
heat transfer rates and heat generation dictate whether a 
reaction will proceed to extinguishment or combustion. 

There are two time constants to be considered in a combustion 
problem of this type: a time constant for the Arrhenius 

reaction and a time constant associated with momentum 
transport. The momentum time constant affects the rate at 
which heat is convected (removed) out of a differential 
volume. The reaction time constant affects the rate at 
which heat is generated (added into) within a differential 
volume . 

D. POROSITY ANALYSIS 

In this problem there are many parameters one might vary 
in order to examine resultant system behavior. Due to time 
limitations this investigation examines the effect of one 
parameter, porosity, on system behavior. The observed 
effect of porosity values ranging from 0.476 to 0.90 is 
discussed. Of the nine runs attempted, only five had suffi- 
cient output data that allowed for comparative analyses. 
Employing Equation 3.1 for the porosity associated with 
spherical particles, the carbon diameter, d, was varied to 
achieve various values of porosity while holding the parti- 
cle center-to-center distance D (Figure 3.1), fixed at 
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0.000417 ft. (0.005 in.). Values of d (ft.) equal to 
0.000417, 0.000411, 0.000396, 0.000385 and 0.000370 were 
used to give values of porosity of 0.476, 0.5, 0.55, 0.60, 
and 0.65, respectively. The results are graphically pre- 
sented in Figures 5.1--5.4. Figures 5.1 and 5.2 show the 
carbon temperature profile and the oxygen concentration 
profiles for porosity p^ = 0.476 (D = 0.000417 ft. = 0.005 
in.). In all cases, excitations for these runs are 1500 

Btu/ft^-hr (in) at — = 0.0 and 50 Btu/ft^-hr (out) at 

o 

— = 1.0; the pressure is 14.65 psi at = 0.0 and 14.55 
r ° z Z ° 

at — = 1.0; the cylinder length to diameter ratio is 0.5 
o 

(length = 1.0 ft.). Figures 5.1 and 5.2 represent typical 
graphical results for carbon temperature and oxygen concen- 
tration profiles. Figures 5.3 and 5.4 are a tabular summary 
of the results of the carbon temperature profile and the 
oxygen concentration profiles, respectively, for varying 
porosities. Profiles with higher values of initial porosity 
are observed to have accelerated development of temperature 
profiles (i.e., as the fuel diameters decrease (increasing 
porosity) , the carbon temperature responds at a faster rate 
to system excitations) . The oxygen concentration response 
is as expected (i.e., as reaction rate goes up, the oxygen 
concentration goes down) . In the runs made for porosity 
values of 0.70, 0.75, 0.85 and 0.90, numerical difficulties 
were encountered. The temperature profiles were observed to 
increase sharply compared to lower values of porosity (300-400° 
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2 

F in the first minute) at the — = 0.0 boundary. Due to 

o 

the steep gradients, further system behavior could not be 
analyzed since output data consisted of a single output 
(T = 5-20 S., problem time, depending on porosity value) 
in 15 minutes of CPU time. 

E. CONCLUSIONS 

The combustion analysis program (CAP) is a viable tool 
for the analysis of heat transfer in porous media. The 
model constructed provides the user with the flexibility to 
solve problems with or without combustion. The role of 
the permeability model must not be underestimated. For it 
is the physical assumptions in this aspect of the model that 
govern the flow, heat transfer and in the end the evolution 
of the combustion problem. 

F. RECOMMENDATIONS 

Follow-on work is recommended in the following areas: 

1. Develop a restart capability for the program which 
will allow for a problem to begin at any point in 
time with initial conditions and rate information 
being identical to previous timestep values. 

2. Nondimensionalize equations in terms of other system 
quantities in addition to spatial dimensions. 

3. Flex the model under other system parameters to examine 
effects on system behavior. 
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4 . 



Consider an internal (particle-to-particle) radiative 
heat transfer analysis. 

5. Consider other fuels and more detailed chemical 
kinetics chains. 
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Figure 5.2 



Oxygen Concentration (Porosity = 0.476) 
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Figure 5.3 Maximum Carbon Temperature Summary Porosity: 
Pi = 0 -476, p 2 = 0.5, p 3 = 0.55, p. = 0.60, 

P 5 = °- 65 - Time: t. = 0.1 minute, 
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Figure 5.4 Oxygen Concentration Profile Summary 

Porosity: p. = .476, p 2 = 0.5, p_ = 0.55, 
p^ = 0.60, p,. = 0.65. Time: t^ = 0.1 
minute, t 2 = 1 minute, t^ = 2 minutes, 
t^ = 4 minutes 
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APPENDIX A 



FORMULATION OF FIELD EQUATIONS 

A. PRESSURE DISTRIBUTION EQUATION 

Darcy's law for two-dimensional flow is. 



Q 



m 3p 
y v 3r 



r + ( 



3p 

3z 



- P 



p 3 g) z) 

d 



( A . 1 ) 



where Q r and Q z are expressed as follows. 



Q 



r 



m _3p 
y 3r 



m _3P 
y 3z 



p p a g) 

cl 



(A. 2) 



(A. 3) 



Invoking the Dupuit-Forcheimer relation, and solving for the 
pore velocity components u (radial velocity) and v (axial 
velocity), Equation A.l becomes. 



-m 3P 
py 3r 



(A. 4) 



v 



-m . 3P 
py 3 z 



P P a g) 

d 



(A. 5) 



The continuity relation (derived in Appendix B) is. 



D(pp a ) 

Dt 



+ p Div V 



0 



(A. 6) 
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Substituting Equation A. 2 and Equation A. 3 into Equation A. 6 
yields , 



3p p 



3t = P p a 5 ' ( ^ ( H r + ( H ' p p a 9)z) 



+ JS.(|£ ; + (|i - pp g ) z ) ' Vp p 
py 3r 3z K a ^ ^a 



(A. 7) 



Expanding terms, Equation A. 7 becomes, 



+ + i-L 

- 2 ~ 2 P 

dr dz a 



+ (• 



3p a 


+ I 


8m 


1 


3y 


3r 


m 


dr 


y 


9r 


9p a 


1 


3m 


l 




Tz~ 


+ — 

m 


3z 


y 


3? 



1. 3P 
r 9r 



3P 

aZ 






p m 

a 



3PP. 

C 

IF 



(A. 8) 



Equation A. 7 with associated boundary conditions (presented 
in Section II. B) is cast into a finite element formulation 
and becomes one of four field equations. Pressure gradient 
information is then substituted into Equation A. 4 and Equation 

A. 5 to obtain the pore velocity components. 

B. POROUS SOLID HEAT TRANSFER EQUATION 

In performing energy balances on both the porous solid 

and on the air, a differential volume of porous medium may 

be partitioned into respective volumes, dV = (l-p)dV for 

s 

the solid, and dV = pdV for the air (shown in Figure A.l) . 

0 . 

The convention used for the energy balance of an arbitrary 
differential volume, dV, is. 
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Heat into DV + Heat generation = Heat out of DV 

+ Increase in internal energy (A. 9) 

The heat transfer mechanisms considered for the carbon 
are conduction, radiation heat transfer between particles, 
convection heat transfer from the particles to the air, 
and heat generation. Applying the above convention, the 
energy balance on a differential volume of porous solid is 
(invoking Tayler Series expansions and neglecting higher 
order terms) , 

[- | f?<(l-p)rq cond _ r ) - fj(U-P)q cond , z ) 

- i(f F (a-p)rq ra d ir ) - !?< < l‘P) q rad , 2 U <JV 

' ^oonv dA ' + = (1 ‘ p), 3int dV (A ' 10) 

In vector form. Equation A. 10 becomes, 

-V* (q + q ,)dV - q dA' + q dA' = (l-p)q. dV 

^cond ^rad ^conv ^gen c unt 

(A. 11) 

Substituting the following expressions into Equation A. 10, 

q , = - k VT Fourier's law (A. 12) 

^cond e c 
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Fourier analogy 
(radiative) 



(A. 13) 



q , = - k VT 

^rad r c 



q = h(T -T ) Newton's Cooling Law (A. 14) 

conv c a. 



'gen 



L int 



= R 



p c C c 3 1 



Heat generation 



Internal energy 



(A. 15) 
(A. 16) 



yields , 



, „ 3T 3T 

i{| 7 (r<l-p)(k e ) + ^(U-PXV 



) }dV 



h ( T - T ) dA ' + R dA ' 
c a g 



3T 



<1_P)P C C c3t dV 



(A. 17) 



Dividing through by dV, and defining dA'/ciV as Z, the specific 
internal area (i.e., surface area per unit volume). Equation 
A. 17 becomes, 



((l-p)(k )VT ) - 

e c 



hZ (T - T ) 
c a 



R Z 

g 



U-pJPcCc 



(A. 18) 



The expressions used to obtain values of the properties and 
parameters in Equation A. 18 are presented in Section III.E. 

C. AIR HEAT TRANSFER EQUATION 

The formulation of the air heat transfer equation begins 
with the general two-dimensional energy balance equation. 



★ 

The difficulty in obtaining an expression for k r is 
addressed in Section III.E. Throughout this work, one may 
keep in mind that k e should really be (k e +k r ) to account 
for conduction as well as radiation within the porous medium. 
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p p D ( U + K) 
£1 

Dt 



= -V*pq - V • pPV - (V'(pi'V)) 



+ hz (T -T ) + p p (V-g) 

C a. 3 



(A. 19) 



where U is the internal energy, K is kinetic energy, x is 

the dissipation function, and g is the gravitational acceleration. 

Expansion of Equation A. 19 yields. 



D , 1 . 2 2 . 

p e aDt <e + 2 (u +v >> " 



-V* ( -pk VT ) - V* (p P V) 

3 . 3 



- I l 

i j i 



, , . p T. . V . + hz (T - T ) 

r f 3x. c 13 3 c a 



+ P P v g 

a. 



(A. 20) 



where 






V(PT-V) = I l 



£ 3x~ p T ij V j “ 7 ^ PT rr V r + pT rz V z 



+ px V + pX V ) 
^ zr r ^ zz z 



| 7 (px rr u +px rz v) + fi-(PT zr u +px zz v) (A. 21) 



And so Equation A. 19 becomes. 



r De D , 1 , 2 2. . 

pp a [ Dt + dX ( 2 <u +v >> 



= V* (pk VT ) -r-(px u+p X V 
^ a a 3r ^ rr ^ 



rz 



_ Tz (pT zr U + pT zz v) " (VpP- V + pPV* V) 



+ hZ (T -T ) + p p v g . (A. 22) 

C 3 3 
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As in the porous solid heat transfer equation, the time and 
position dependent porosity appears inside the differential. 
Expanding terms, the air energy equation becomes, 

p P a§t + ~T' §t (u +V ) = (P k a VT a ^ “ [V p P ' V +p PV-V] 

g 

- -r— (pT U + pi V) 

9r ^ rr ^ rz 

- l—(px u+px v) + hZ (T -T ) 

3z ^ rz ^ zz c a 

+ p p v g . (A. 23) 

3 . 

Consider the momentum equation for the r-direction. 



R Momentum 



3t (p p a u) 



JL / Om £ % Of i 

F ( - 3X (r P p a u > - Tz P p a uv 



3pP 

9r 



,13 , . PT 

~ ( ? 9? (p rT rr ) ~ — r 



0 0,3 , . . 

+ 37 (P T rz 



(A. 24) 



Consider the momentum equation for the z-direction. 



Z Momentum 



3t ( P p a V) 



F ( - lx (r P p a uv) - !x ( P p a v2) - ll ( P P) 



- ( F f? (p T rz> + Ix'P T zz n + p p a 9 



(A. 25) 
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and the continuity equation 



Continuity 



3t'P p a> 



= - u 35 ( P p a> 



(P P 



3u 
a 3r 



+ 



P p a U , 



V 



3pp a 

3z 



- P P 



a 



3v 

3z 



(A. 26) 



Multiplying the continuity equation through by u, and substi- 
tuting this into Equation A. 24, the r-momentum equation 
becomes. 



P P 



3u 
a 3t 



3u 

pp a u 3? 



P P v 

U 



3u 

3z 



3pP 

3r 



.13 , . 

- ( ? S (pn rr ) - 



px 



69 



9 t 

37 (p 



rz 



) ) 



(A. 27) 



Multiplying Equation A. 26 by v, and substituting this into 
Equation A. 25, the z-momentum equation becomes. 



3v 

pp a 3t = 



3v 
'a u 3r 



3v _ 3pP _ ( 1 3 ( 
a 3z 3z 



- P P= - p p, v — - - (- ^(prx^) 



3 , . 

+ XT P x ) ) + p p g 
9z zz a 



(A. 28) 



Multiplying Equation A. 27 by u and noting that. 



Du 

P p a u Dt 



3u 



2 3u 



3u 



= PP a u H +Pp a U 37 + p p a UV 37 



(A. 29) 
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Equation A. 21 becomes 



Du 

P p a U = 



3pP 
3 r 



/I 3 / 

’ U ’ u( r 3T (P rT 



rr 



PX 90 , 3 , 

— + I? (P T rz 



(A.30: 



Multiplying Equation A. 28 by v and noting that. 



Dv 

P p a V Dt 



3v 3 V 2 3v 

P p a V 3t + P p a U V 37 + P p a V Ti 



(A. 31) 



Equation A. 28 becomes 



Dv 

P p a V Dt = 



- v 



3pP ,13, , 

—— - v (— tt— (p rx ) 
3z r 3r ^ rz 



+ 3l ( P T ZZ >> + P P a V 9 



(A. 32) 



Expanding Equation A.30 yields. 



P p a U 



Du 

Dt 



3pP 9pT rr 

u - u — r 

3r 3r 



u x + U P T 9 8 
r P T rr r 



u 3l ( P X rz ) 



(A. 33) 



Expanding Equation A. 32 yields. 



Dv 

P p a V Dt = 



3pP 3(pT rz' 

' v 3z " V 3r 



FP T rz ' V h ipT zz ] + ? P a V g 



(A. 34; 



The energy Equation A. 23, after substituting and expanding 
terms, is. 
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De . Du , Dv 

P p a Dt + p p a U Dt + p p a v Dt 



= V -(pk a VT a ) - (P T rr |f 



3<pT rr' , 3v 

+ u —3? + pT rz 3? 



8 , x . , 3u 

+ v 3? ( P T rz>> ■ P T zr 3? 



3 , . , 3v 

+ u P T :r + pT zz 37 



+ V -r— (px ) ) - ( VpP • V 

dZ ZZ 



+ p PV-V) + p p V g 

3. 



+ hZ (T - T ) 
c a 



(A. 35) 



Substituting Equation A. 33 and Equation A. 34 into Equation 
A. 35 yields. 



Dp -> 

p p r— = V* (pk VT ) - p P Div V + hZ(T -T ) 

3 Ut 3 3 C 3 



, 3u 

(p T rr'S7 + 



P T rz37 + 



3v 

TT + P ^ + 



zr 



3u 

3z 



3v. 

^ T zz 3z 



u p x 



99 



u 

+ — n T 
r F rr 



+ v pT 

r F rz 



(A. 36) 



The viscous dissipation terms in Expression A. 36 are 
neglected because the fluid is a gas flowing at low velocity 
(see [Ref. 9]). Therefore, the energy equation for the air 
in the porous medium is. 



95 



(A. 37) 



P P 



De 
a Dt 



V* (pk VT ) - pPDiv V + hZ(T -T ] 

3 3 C 3. 



With specific enthalpy for a gas defined by. 



h = e + - 

P 



(A. 38) 



De/Dt can be expressed as, 



De 

Dt 



Dk 

Dt 



1 DP 
p a Dt 



Dp. 



2 Dt 



(A. 39) 



where k is the specific enthalpy, and e is the specific 

energy (internal plus kinetic). Multiplying the continuity 

2 

equation through by P/(pp ) obtains, 

3 



P 3 



P 1 3 



PP 



2 Jt (pp a> + — 7 3? (p p a ru) + ~ 3l lpp a v) 
a PP = 



_P_ 9 

pp; 



= 0 . (A. 40) 



Expanding yields. 



PP. 



, ^ p a 3p. 

(p — + p a 3t> 



P 3 , 

■2 u 3? <pp i 



pp a 3 , , 

T-(ru; 



PP: 



PP. 



2 r 3r 



P v 3_ 

2 3z 
PP = 



( P P a ) - 



P 3 v 

2 P p a37 
PP = 



(A. 41) 



In vector notation, Equation A. 41 becomes. 
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(A. 42) 



_ + _p_ i£ = 

2 9 1 PP 3t 

p a a 



4(V-7) p - -£-(V*7)P 

2 pp_ 



— Div V 
p a 



Employing Stokes (substantial) derivative notation. Equation 
A. 42 becomes. 



Dp 

P H a 
Dt 
p a 



-L- S - 1 Div V 
P p a Dt p a 



(A. 43) 



Substituting Equation A. 43 into Equation A. 39 yields. 



De _ Dh 1 DP 
Dt Dt ~ p Dt 

a. 



Div V 



P Dp 
' pp a Dt 



(A. 44) 



Substituting Equation A. 44 into Equation A. 37 yields. 



Lk DP _ Dp 

DO " P “ P — 



Dt Dt 



V* (p k VT ) + hz (T - T ) 

a a c a 



(A. 45) 



Simplification yields, 



P P. 



Dh _ DpP 



aDt Dt 



V* (p k VT ) + hZ (T - T ) 

a a c a 



(A. 46) 



Invoking the Maxwell relations for a simple compressible 
substance , 



dh = T ds + - dP 



(A. 47) 
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ds 



(A. 48) 



= c 4? - 1 dP 

P T p 



and recalling that. 



6 = - 1 !£ ) 

^ p 9T p 



the equation of state for a perfect gas. 



P = p RT 



simplifies Equation A. 49 such that. 



B = - ££(-=%) - i 

RT 



Thus, Equation A. 48 becomes. 



ds = C ^ 7p dP 

p T pT 



Substituting Equation A. 52 into Equation A. 47 and 
terms yields, 



• dh = 



C dT 
P 



or 



Dk 

Dt 



= C 



DT 

a 

a Dt 



(A. 49) 

(A. 50) 

(A. 51) 

(A. 52) 
cancelling 

(A. 53) 
(A. 54) 
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Substituting Equation A. 54 into Equation A. 45 yields. 



3T 



PP = C 



a a at 



V • (pk VT ) - pp C (V-V)T 

3 3 3 3 3 



+ hz<T c -T a ) + 



(A. 55) 



Making the assumption that pP changes very little with time 
[Ref. 9], (this assumption was subsequently confirmed by the 
model) , i . e. , 



DpP 

Dt 



(V- V) pP 



S <P p ) + v 3 <P p > 

or dz 



(A. 56) 



The final air heat transfer (energy conservation) equation 
is. 



3T 



PP C 



a a 9t 



V* (pk VT ) - p p C (V-V)T + u + v 

^ a a ^ H a a a dr oz 



+ hZ (T - T ) 
C 3. 



(A. 57) 



In vector form. Equation A. 57 becomes. 



9T 



PP C 



a a at 



V* (pk VT ) - p p C (V* V) T 

3 3 3 3 3 



+ (V-V)pP + hZ(T - T ) 

C 3 



(A. 58) 



The expressions used to obtain the properties and parameters 
in the coefficients of Equation A. 57 are presented in 
Section III.E. 
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D. OXYGEN MOLECULE DIFFUSION EQUATION 

The final consideration in formulating the field equations 
for the model is the transport of oxygen molecules. The 
oxygen molecule transport equation is obtained by a conser- 
vation of species balance on the differential volume of air, 
dV = pdV. The convention used for the species balance into 
a differential volume is, 



O 2 in = C >2 out + 0 2 consumed + 0 2 accumulated (A. 59) 

The transport mechanisms considered were diffusion due to 
concentration gradients (Fick's law), convection, and air 
consumption by combustion. The species balance on oxygen 
becomes , 



p”diff dA 



+ pm diff dA 



+ pm dA 4- pm dA 
c conv c conv 

z r 



* pm diff dA 



r+dr 



+ pm diff dA 



z+dz 



+ pm dA 
conv 



r+dr 



+ pm dA 
conv 



z+dz 



+ m dA' + pm dV 
cons acc 



(A. 60) 



Representing terms on the right side by Taylor series expan- 
sions (neglecting higher order terms), i.e., 



pm.dA = pm.dA + - (pm,dA)de. 

K . j K <3 £ • K 1 



£ i +de i 



(A. 61) 
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where , 



k = conv, = r or z dA r = r d0 dZ 



(A. 62) 



and dA z = r dr d0 



(A. 63) 



Then Equation A. 61 becomes, 



■f^:(p mdA r )dr = ^(r pu <J>d0 dz) dr 



p |^(r P u <J>) r d0 dZ dr 



(A. 64) 



and 



■^■(p mdA z )dz = -^-(p v <J> rdr d0)dz 



= j^(p vi|))r d0dz dr 



(A. 65) 



Substituting Equations A. 64 and A. 65 into Equation A. 60, 
cancelling terms and rearranging. Equation A. 61 becomes. 



3 * 3 * 

-r— (p (m, . + m ) dA ) dr - -r— (p (m, . + m ) dA ) dz 

3r ^ diff conv r 3z ^ diff conv z 



- m dA ' = pm dV 

cons acc 



(A. 66) 



Substituting the following expressions into Equation A. 66, 
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(A. 67) 



m , . ^ ^ = - V V <p 

diff e v 



m = utb. m = vd> 

conv,r conv, z 



(A. 68) 



m = R~ 

cons 0, 



(A. 69) 



m 



acc 3t 



(A. 70) 



dividing both sides by dV, and letting dA ' /dV equal the 
specific internal area, z, the oxygen molecule diffusion 
equation becomes, 

V- (pP e Vcf>) - V- (p V) - R Z = (A. 71) 

The methods and expressions for obtaining the properties and 
parameters in the coefficients of Equation A. 71 are presented 
in Section III.G. 
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POROUS MEDIUM 



AIR 



SOL 1 0 




Figure A.l 



Separating A Differential Volume of 
Porous Medium into Respective Volumes 
of Solid and Air 



APPENDIX B 



AUXILIARY EQUATION FORMULATION 
A. CONTINUITY EQUATION 

The continuity equation for a fluid in a porous medium 
is expressed by, 

D(PPJ 

— p ~— - + P P a Div V = 0 (B . 1 ) 

or in equivalent form, 

3 

■g^(pP a ) + (V- V)pp a + pp a (V*V) = 0 (B . 2 ) 



and 



3 

3t 



(pp ) + u 

d 




V 



9 (pp_) 

d 

9"z 



,1 

pp a ( ? 



9 (ru) 
3r 




(B . 3 ) 



B. LAGRANGE POLYNOMIAL APPROXIMATIONS FOR THERMAL PROPERTIES 
Relations for calculating the dynamic viscosity, thermal 
conductivity, and specific heat at constant pressure of air 
at different temperatures were required. Second order 
Lagrange polynomial fits to empirical data provide a simple 
method for obtaining the relations required. Vatikiotis 
[Ref. 9] gives the details for arriving at the resulting set 
of polynomials, 
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y = -3.308 xio" 9 T 2 + 4 .633 x 10~ 5 T + 4.427 x 10 2 (B.12) 

3 . 



k a = -2.608 xio' 10 T 2 + 1.930 x 10~ 5 T a + 1.361 x 10 2 (B.13) 

C = -1.293 xio' 9 T 2 + 2.758 xio" 5 T + .238 (B.14) 

3 3 3 



Each expression obtains property values within two percent 
of the data presented in [Ref. 9] for temperatures up to 
3000 degrees Fahrenheit. 

C. CARBON PARTICLE SURFACE RECESSION 

The following analysis of particle diameter consumption 
assumes that the fuel particle surface recedes uniformly. 

This assumption is reasonable if the particle diameter is 
small in comparison to the cylinder geometry, i.e., negligi- 
ble boundary effects. In addition, since the velocities are 
low, the hydrodynamic effects on the uniformity of the 
particle surface recession are negligible. The analysis also 
assumes there are no significant thermal gradients within 
the particle, i.e., an isothermal carbon particle. This is 
also reasonable for the small particles examined in this 
analysis. A mass balance must be performed on the particle 
and equated with the reaction rate. This equivalence yields 
the following expression, 



= R Z 



c 



dm 

dt 
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(B . 15 ) 



or equivalently. 



d (1-p) p 

c 

dt 



R Z 
c 



(B. 15a) 



Simplification of Equation B.15a yields, 



dp 

dt 



V 



(B. 16) 



Substituting the porosity expression for spherical particles 
yields , 



— (-(-) 3 ) 
dt '6 'D' ’ 



R Z 
c 



(B . 17 ) 



Application of the time operator yields. 





R Z 
c 




(B. 18) 



Isolating the diameter time derivative yields, 

- 2 R ZD 3 

d = (B . 19) 

P c TTd 
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APPENDIX C 



GALERKIN FEM FORMULATION 
A. FINITE ELEMENT METHOD 

The solution of the system of coupled, nonlinear partial 
differential equations given by Equations 3.9, 3.22, 3.26 
and 3.27 — subject to boundary and initial conditions, was 
obtained by a Galerkin formulation of the Finite Element 
Method. 

1 . Galerkin Formulation 

A Galerkin formulation of the Finite Element Method 
was used to obtain solutions of the porous solid and air 
energy equations, the oxygen diffusion equation, and the 
continuity (pressure — Darcy's law) equation. A convenient 
form of Equations 3.9, 3.22, 3.26 and 3.27 was used in the 
formulation where the spacial coordinates, r and z were 

nondimensionalized by e = z/z and n = r/r . 

2 o o 

The closed domain defined in (r,z) space by (0,0), 
(0,1), (1,1), and (1,0) was partitioned into NEL 

(2* (NRNP-1) * (NZNP-1) ) contiguous area elements obtained 
from a NSNP model. NSNP , NRNP and NZNP are the number of 
system nodal points, radial points, and axial points, 
respectively. The four field variables T , T , P and 
were approximated by, 
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NSNP 

T = iK.(n,£,t) = l N . (n / £ ) 6 . . ( t) (C.l) 

c - L J i=l -*■-> 



T a = ^ 2 j (n ' £ ' t) 



NSNP 

l Nj (n»e) 9 2 j 



( C . 2 ) 



NSNP 

P = <K.(n,£,t) = l N. (n,e) e, . (t) (C . 3) 

J 3 j=l 3 J 3 



<P = <l< 4 j (h /£ ,t) 



NSNP 

I N ( n/£ )e (t) 

j=l 3 4D 



( C . 4 ) 



where N^ for j = 1,...,NSNP is a set of specified linear basis 

functions with local support, and the sets 9^, ®2j' ®3j' 

and 9^^; j = 1,...,NSNP, are the solution coefficients to be 

determined. The N. were selected to satisfy the condition 

N.(NP.) =6.- where the Kronecker delta, <5.., is defined by 
3 i 13 13 * 

6 . . = 1 for i = j , and 6 . . = 0 for i 4 j . As a result, 0. . , 

i-j j j lj 

© 2 j / 9^j and 9^ are the values ip ^ , \p 2 , \p ^ and \p ^ at the nodal 
points (i.e., i[Kj(ri,£,t) = 9^j (t) ) . 

Area interpolation functions (shown in Figure C.l) 
were used as the linear basis functions which provide the 
necessary function continuity. As a measure of error, a 
residual, , is defined for each field equation by. 



R. 

l 



= ip + A. (ip) - F + C [T <f> ] * 
~ 1 



(C . 5) 
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where denotes the spatial operator of the ith equation 
and the asterisk denotes that this term appears in the carbon 
temperature and oxygen concentration equations only. The 
term arises from the reaction terms and is developed in 
subsection 3 of this section. The following notational 
convention for differentiation is adopted, 



For field equations 3.9, 3.22, 3.26 and 3.27, the residuals 
are , 




H ) 

3i 



(C . 6 ) 



3 ( ) 

at 



(C . 7 ) 



Rt = ( 1 -P)P c C c T c - V- ((1-p) (k e )7T c ) + hZ(T c -T a ) 



R Z 

g 



(C . 8 ) 



R, 



T 



A 




+ P P a C a (V ’ V)T a 



(V* V)pP 



(C . 9 ) 




(C. 10) 



R 



<J> 



p <f> - V ' (p V V <p) + V • (p <j> V) + R n Z 

S U2 



(C.ll) 
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where the coefficients of the response variables are them- 
selves functions of the response variables, and thus, the 
equations are nonlinear. In accordance with the Galerkin 
method, the final system of ordinary differential equations 
was obtained by setting each residual, R^, orthogonal to 

each basis function, N., that is, 

3 



/ / N R dA = 0 (C . 12 ) 

A " 

e 

The 4*NSNP ordinary differential equations given by Equations 
C.12 retain the character of the original set of partial 
differential equations, i.e., self-adjoint operators yield 
symmetric matrices and non self-adjoint operators yield 
nonsymmetric matrices. Thus, linear field operators trans- 
form to matrix operators and nonlinear, coupled algebraic 
operators. Incorporation of the boundary conditions resulted 
in 4*NSNP nonlinear coupled ordinary differential equations, 

F(ip,ip,t) = A \jt + Bf - F + [i^ T ip^] (C . 13 ) 

subject to initial conditions, where A is a (4*NSNP) * (4 *NSNP) 

matrix, B is the matrix associated with the linear field 

operators in Expression C.5, F is an excitation vector, 

* 

is a 3x9 matrix arising from a bilinear operator treat- 
ment of the reaction terms in Expressions C.8 and C.ll. 
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Adopting the convention, 




NSNP 

I N.0 . . , i = 1,4 

j = l ^ ^ 



( C . 14 ) 



and performing an integration by parts on the second order 
derivatives yields, 



//NR dA=-£ (l-p)N(k ) (T ) d£ + // 

A c Ze A 

e e 



(1-p) (k ) N (T ) dA 
^ e ~r c r 



- // N(k )dA + // 

A A 



(1-p) (k ) N (T ) dA 
e ~z c z 






le 



(1-p) (k ) N ( T ) dA + / / hZN (T -T ) dA 
e-cn J J ~c a 

e 



- If 

A 



R ZNdA + // p p NT dA 



(C . 15) 
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^ = ^P k a? (T a ) n d£ + // pk a5?r (T a } r dA 
3- .A 



pk 

- // — - N(T ) dA + // pk N (T ) dA 
{ J r ~ a r { J ^ a- z az 

A A 

e e 



-6 pk N(T ) dS, - // hZN ( T -T ) dA 
£e a ~ a n A ~ c a 

e 



- ff upN(P) r dA - 



// u NPdA 

A 9r " 
e 



- // vpN (P) dA 
A ~ z 



- // 

A 

e 



v 



NPdA 

a Z ~ 



+ // P P a C a uN(T a ) r dA + ff p p a C a vN(T a ) z dA 



// 

A 

e 



p C N 
*a a~ 



dA 



(C. 16) 



ffl 



N RdA 
P 



p m 

-4 -^N(P) dZ + f f -2— N (P) dA 

n A y ~ r r 



p a m 



- ff — N (P) dA - 6 P * m M / D x 
V yr ~ r r n „ — N(P) d£ 



Ze y n 



+ // 



P,m m 3pp a 



} J —— N (P) dA + // — *-=• N (P) dA 

A y ~z z a ; py °r ~ r 



3pp 

+ // & S (p) z dA + // sf- 5 pdA 

A. A 3. 

e e 



(c. 17; 



//NRdA = -^ te pP e N(d)) n dA + // P P e N r (4>) r dA 



pP 

// 6 



N ( 4> ) dA 
r ~ r 



4 pP N(4>) d£ + // pP N (<f>) dA 
e ~ n A e ~ z 2 



(rpu) 



+ // puN(<J>) r dA + // pvN(4>) z dA + // - — - N<J>dA 



+ // (pv) N<f>dA + // R ZNdA + // pN<f>dA 

A o 2 ~ A o ~ 

e e e 



(C . 18 ) 



The line integral terms in each expression above are boundary 
terms which permit incorporation of natural boundary conditions. 
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Implementation of the boundary conditions is presented in 
Section 2 of this appendix. The coefficients in Equations 
C.15, C.16, C.17 and C.18 are temperature dependent properties, 
and were taken as the average values of the properties over 
an element. In the limit, as the elements get smaller (i.e., 
NSNP -+■ °°) , the average values of the coefficients converge 
to the exact values. Inspection of expressions C.15, C.16, 

C.17 and C.18 yields the six operators, 




(C . 19 ) 




(C. 20) 



// N OjO dA 
A r 

e 



(C. 21) 



If N „(i|0 dA 



-z z 



(C.22) 



// NijdA and // Nt^dA 

A ~ A ~ 



(C. 23) 



N (ip) dZ 
Ze~ n 



(C . 24 ) 
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To formulate these operators, the global linear 



shape function, , are defined on the local level by, 



T T 

N 0 -+ £ 6 



[C. 25) 



where the natural coordinates ^ an< ^ £3 are defined by 

Figure C.2, 



A. 

q = ^ i = 1/2,3 (C . 26 ) 

e 

and A is the area of the elemnet and the A. are the areas 
e 1 

(in Figure C.2) subtended by lines from a point P(r,z) 
inside the triangle to the triangle's vertices (rj,Zj, 
j = 1,2,3). The local shape functions (i.e., the elements) 
have the following properties. 



^i j = 0 — • ; i = 1,4, j = 1,3 (C. 27) 

£ .(NP.) = 6 . . ; i = 1,3, j = 1,3 (C.28) 

-L J ± J 



Having defined the local shape functions, the ele- 
mental matrix operations C.19 through C.24 are, 



// N (ifO r dA 6 1 ®elt 

A 

e 



// N(») dA - |lc j ]9 elt 
A 

e 
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(C. 30) 



1 



(C . 31 ) 



// ^) r dA 



e 



4A 



■^iVSeit 



// N (ip) dA = 



A 



4A 



• [c i c j ] ®elt 



A 

jj N^dA = j|[k]0 

A 



elt ' 



k = 1, i jt j 
k = 2, i = j 



6 kN (\p) dl = kN(tJj) dr - <j> kN(^) dz 

Jle ~ n £e ~ r £e ~ 

where a(i,j) coefficients of the element matrices 
given by, 



a ( i , j ) = [g] ; g = g(i, j) 



and the vector is given by, 



.elt 




The derivations of these operators are presented in 
A. 5 of this appendix. 



(C . 32 ) 



(C. 33) 



(C. 34) 



(3x3) are 



(C . 35 ) 



(C. 36) 



Section 
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2 . Implementation of Boundary Conditions 



Having formulated the system matrices for the field 
equations, treatment of the boundary conditions is now dis- 
cussed. Each field equation is treated individually. 

a. Porous Solid Energy (Heat Transfer) Equation 

There are three types of boundary conditions that 
can occur: Dirichlet, Neumann and Cauchy (mixed) boundary 

conditions. The treatment of each is as follows. For 
Dirichlet boundary conditions, one may specify an equation 
to be a linear equation (i.e., independent of time) by placing 
a 1 on the diagonal of the system stiffness matrix correspond- 
ing to the particular degree of freedom at hand and setting 
each time derivative matrix coefficient for the same equation 
equal to zero. An alternative scheme involves setting the 
time derivative diagonal term equal to one and setting all 
stiffness coefficients (in the same equation) equal to zero. 

In this manner one specifies the Dirichlet value as an initial 
condition whose residual is identically zero and thus is 
invariant with time. Treatment of Neumann boundary conditions 
is as follows, 

4 kVxd£ = 4 pdf; -»■ — r— added to local F(l) 

Ze He 

and local F(2) at the element 
nodal points (C.37) 

In the above expression, x i- s a response variable and p is 
an average value of flux across an element. Treatment of 
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Cauchy boundary conditions is as follows, 



j> kVydC 
£e 



•f> Mx -xj ■* x le -r 

le 



2 1 
1 2 



X 1 

X 

2 i 



— 0 

scattered into 8 and - K, — added to 



F(l) and F(2) at the element nodal points 



(C. 38) 



For a complete development of the theory for incorporating 
boundary conditions, see [Ref. 28]. 

3 . Treatment of the Reaction Rate Term 

An exponential reaction rate term appears in both 
the porous solid and oxygen diffusion equation. The reaction 
expression is, 



R c = A <J) n exp(-E/R u T c ) 



(C . 39 ) 



In the carbon equation, Expression C.39 appears as. 



r Z = (R AH )Z 
a r. K 



(C . 4 0 ) 



In the oxygen concentration equation, Expression C.39 
appears as. 




<Vr 1,z 



(C . 4 1 ) 
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The reaction rate term may be treated in various ways, 

a. Treatment as an Excitation (Force Term) 

This treatment is the easiest of all schemes. 

The term is merely evaluated at the last time step and is 
assumed to be constant over the next integration time step, 
i . e . , 



A ^ 

R = {A exp (E/R T ) } (C.42) 

G U2 LI C 

The superscript * indicates evaluation occurring at the 
previous time step. The value of the term is evaluated 
at the ith nodal point and inserted as a factor in the 
corresponding carbon energy or the oxygen diffusion equation 
of the ith nodal point. 

b. Linear Operator Treatment of O 2 Concentration 

In order to realize an improvement over the first 
treatment, one may retain a portion of the reaction expression 
as an operator by making the following rearrangement: 

R = (A R^ <{> n ~ 1 exp(-E/R i T ) }*• [<J>] (C.43) 

G (J ^ LI C 

where [ 4> ] denotes a spatial operator treatment of the response 
variable . 

c. Temperature and 0- Concentration Bilinear 
Operator Treatment of the Reaction Rate Term 

The bilinear operator treatment of the reaction 

rate is the present method of treatment. The reaction rate 
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term is rearranged as follows, 



A R 



n-1 



exp(-E/R u T c ) * 



R 



‘V 1 



(C.44) 



Letting 



T 

c 




and 



<P = 




(C.45) 



(C. 46) 



and invoking natural coordinates [Ref. 28] , an elementally 
averaged contribution results in the following area integral. 



C 



c // £ 

Z A 



e 



T 

§ e 4 dA 



Expanding yields. 



= C II f < 5 1 5 j 5 C 2 §-i UjSj 



li °4 j 



i = 1,3, j = 1,3 



(C . 4 7) 



(C. 48) 



A final explicit expansion of the integral looks like, 
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c = c 



ih 



h4 

h^3 



4h 

hhh 

hi 



• h^2 

•c c 2 

. *1*2 

I ^1^2^3 



C c 2 

n^2 



«2«3 



' - 2 , 
Cl 3')^>'5 * C 1 S 



1*2*3 J *1*3 

: w 3 



*2^3 



c 2 c 3 



2 * r r 2 

i *1*3 



*1*2*3 "1*3 



«2«3 



*2*3 



. r 2 
’2 '”3 
.2 



[C.49) 



Invoking the integral formula of natural coordinates [Ref. 
28] , 



// C* C2 £3 dA 

A 

e 



aiblcJ 2A 



(a + b + c + 2) ! 



(C. 50) 



results in. 
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(C . 51 ) 
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2 


t 

* 


1 


2 


2 


1 

, 2 


2 


6 







where , 



A Z r£ <], n 1 exp(-E/R T ) 

C = < (AH r or f R ~)* 



-1. 1 



(C. 52) 



The asterisk denotes as applicable for T c or <j> equations. 

4 . Implementation of Reaction Term in the Numerical 
Me thod 

Franke's (modified Gear) integration routine requires 
a calculation of (or an approximation to) the Jacobian matrix 



121 



from the system matrices, A(t), B(t) and C(t). In order 
to include effects of T c and t}> arising from reaction rate 
terms in the Jacobian vector, and thereby improve the effi- 
ciency of the integration routine, the combustion terms are 
incorporated into the residual equations (in DIFMOD) . 
Modifications to the Jacobian matrix are accomplished in the 
JACMOD and NUITSL routines. Reaction rate terms of Expres- 
sions C.40 and C.41 contribute nine terms to the respective 
residual equations and twenty-seven terms to the Jacobian 
vector. Reaction terms are generally expressed as, 

3 3 

RT = l l < “ i ]<; ( ^ li ^4i ) ' k = 1,9 (C.53) 

i=l j=l J 



The Jacobian is defined as. 



3 ^ ^ 



(C . 54 ) 



Contributions arising from RT as a result of dF/dip are, 



3RT 

3 *li 



3 

I 

j=l 



^ik* 




(C . 55 ) 



where the k* are compatible with the column locations in the 
C matrix of the products, and 
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3RT 



3 




i=i 




(C. 56) 



where the k* are compatible with the column locations in the 



may be incorporated into the PW Jacobian vector (containing 
contributions to the Jacobian from the linear spatial 
operators) as contributions from Equation C.8 with combustion. 
Similarly, terms from Equation C.56 may be incorporated into 
PW for terms arising from the C >2 residual equations with 
combustion. The remaining contributions to the Jacobian 
vector, Equation C.56 for the carbon equations and terms 
arising from Equation C.55 for the equations are stored 
in the PWMOD (NDOF/2,9) matrix. The column number or variable 
of differentiation array, INMOD (NDOF/2,9) and PWMOD array 
are communicated to the NUITSL routine so additional Jacobian 
terms not assimilated into the PW array (by virtue of the 
storage scheme selected) may be taken into account during 
the convergence sequence for the new iterates. The iterating 
scheme is of Newton-Raphson type. Thus, the final synthesis 
of the Jacobian including effects of all terms arising from 
combustion is consummated. 

5 . Derivation of the FEM Operators 



the following six differential operators were identified, 




In the section on the finite element formulation 
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(C. 57) 



II N ( ip) dA 
A ^ r 
A e 



II N(4») dA 
A - z 

e 



(C . 58 ) 



II N (<P ) dA 
A " r r 
e 



(C . 59 ) 



// N (ip) dA 
A ~ z z 



(C. 60) 



// N(^) dA and ff N \p dA 



:c. 61 ) 



j> kN (\|) ) d«, 
ie ~ n 



(C . 62 ) 



where N ^ are the global basis functions. These operators 
are constructed on the element level by introducing the 
corresponding element basis functions, ^ . The global and 
element basis functions are related by. 



T T 

= n e -*■ n 



(C. 63) 



The derivation of the local elemental matrices (using 
the local coordinate system depicted in Figure C.3) according 
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to the Galerkin method for the global operations proceeds 
as follows: 

For Operator C.29(also C.57), 

Global Local 



// N(^) dA 

A ~ 
e 



// ? re dA 

A ~ ~ r ~ 



(C. 64) 



Noting that, 



lyr 



b .6 . 
2A 



and 



j 

3 z 



c .9 . 
2 2 
2A 



j = 1, 3 



(C . 6 5 ) 



where the repeated index implies Einsteinian notation, the 
elemental matrix becomes, 



// § 

A 



b T 



2A 



-6dA 



(C.66: 



Expanding , 



// 



2 
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b . T 



e 



6dA = =r 



b, b. 



b, b. 



b l b 2 




6 (C. 67) 



For Operator C.30 (also C.58), 



Global 



Local 



// N (ip ) dA 
A z 



// § §2~ dA 



(C. 68) 
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Substituting the local shape functions gives 



// 



A 

e 




0 dA 



the elemental matrix becomes, 




For Operator C. 31 (also C.59) 



Global 



Local 



// N (i|0 
A ~ r 



e 



r 



dA 

e 




Substituting the local shape functions gives, 



// 



A 



2A 



b. 

2A 



0 dA 



If 



4A e A 
e e 



b l b 2 

b l b 3 



b l b 2 



b 3 b l 



b l b 3 

b 2 b 3 



0dA 



the elemental matrix becomes, 



II N <^) _ dA 

A ~r r 

e 



4A 



b l b 2 

b l b 3 



b l b 2 



b 2 b 3 



b l b 3 

b 2 b 3 



(C. 69) 



(C. 70) 



(C. 71) 



(C . 72) 



( C . 7 3 ) 
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For Operator C.32 (also C.60), 



Global 



Local 



// N U) dA 



z z 



// C z € z 6 dA 

A - z ~ z ~ 
e 



Substituting the local shape functions gives, 



T 

C . C 

[} 2A 2A 
A e e 

0 ~ ~ ~ ~ ~ ~ 



0 dA = 



4A 



If 

A 



e e 



C 1 C 2 

C 1 C 3 



c c c c 

u l 2 ^1 3 



C 2 C 3 



e dA i 



C 2 C 3 



the elemental matrix becomes, 



If c z c z e dA 

A ~ Z ' Z 



4A 



C 1 C 2 

C 1 C 3 



C 1 C 2 



C 2 C 3 



C 1 C 3 

C 2 C 3 



For Operator C.33 (also C.61) 



Global 



Local 



If N ip dA 
A ~ 



II n e dA 

A ~ ~ 



Substituting the local shape functions gives, 



C . 74 ) 



C. 75) 



C. 76) 



C .11) 
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// 

A 



C r 0 dA = // 



^2 

^1^3 



^ 1^2 

^ 2^3 



^ 1^3 

^ 2^3 



a 



0 dA 



(C . 78 ) 



the elemental matrix becomes. 



// 

A 



C £ 9 dA 




(C. 79) 



The last operator C.34 (also C.62) is incorporated into the 
excitation vector as described in the FEM formulation and 
has been addressed in Subsection 2, Implementation of 
Boundary Conditions. 
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Figure C.l 



Area Interpolation Functions 
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Figure C.2 Area Coordinates 
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Figure C.3 Local Coordinates 
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